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Abstract 

We use a numerically exact real-time path integral Monte Carlo scheme to com- 
pute electron transfer dynamics between two redox sites within a spin-boson ap- 
proach. The case of asymmetric reactions is studied in detail in the least under- 
stood crossover region between nonadiabatic and adiabatic electron transfer. At 
intermediate-to-high temperature, we find good agreement with standard Marcus 
theory, provided dynamical recrossing effects are captured. The agreement with our 
£f) ' data is practically perfect when temperature renormalization is allowed. At low tem- 

perature we find peculiar electron transfer kinetics in strongly asymmetric systems, 
characterized by rapid transient dynamics and backflow to the donor. 
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1 Introduction 



This paper presents results of computer simulations for simple asymmetric 
electron transfer (ET) reactions in polar solvents. Such reactions are encoun- 
tered in a large variety of different systems, e.g. charge transfer in semicon- 
ductors, chemical reactions, or the primary step in bacterial photosynthesis 
[1,2,3,4,5,6]. Our theoretical study of such processes is based on the spin-boson 
(SB) model (dissipative two-level system) [7]. This model provides an accurate 
description of many ET reactions involving condensed-phase environments [4], 
if a suitable spectral density J{uj) can be determined for these "bath" modes. 
For a detailed discussion of restrictions for the SB description of ET reactions, 
including an (almost) exhaustive list of references to work done in that context, 



Email address: lothar.muehlbacher@physik.uni-freiburg.de 

(L. Miihlbacher). 



Preprint submitted to Elsevier Science 



2 February 2008 



we refer to our previous work [8] for symmetric ET reactions. In the present 
paper, primary focus is laid on the asymmetric case. Within the SB model, 
the localized sites representing the donor and the acceptor are described in 
terms of an effective, quantum mechanical "spin" variable. The two-level sys- 
tem (TLS) is characterized by a tunnel splitting hA, which is twice the usual 
electronic coupling between the redox sites, and by the asymmetry he giving 
the energy difference between the two localized states. The frequency-resolved 
coupling strength of the electron to the surrounding solvent modes is then 
described by the spectral density, where one mimics the environment by an 
infinite set of effective harmonic oscillators. The spectral density for a given 
ET reaction can in principle be computed from classical molecular dynamics 
simulations [4], but in many cases it turns out to be appropriate to consider 
the special class of Ohmic spectral densities [5,7], 

J(u) = 2itau e~ u ' u % (1) 



with dimensionless damping strength a and a cutoff frequency u c . On general 
grounds, the linear low-frequency behavior of J(co) is expected in basically all 
condensed-phase ET reactions [4,5], and the frequency u c then corresponds to 
some dominant bath mode. To make contact with common notation, we shall 
use the reorganization energy hA [1] instead of a, where A = 2au c . While 
our method below is able to treat arbitrary spectral densities, for clarity we 
will only present data using Eq. (1). Note that for common ET reactions, 
one has large damping parameters a, and therefore one is always in the in- 
coherent regime with respect to the electronic degree of freedom. Coherence 
is nevertheless possible with respect to the bath degrees of freedom, in par- 
ticular when fast bath modes are absent or particular initial preparations are 
selected [9,10]. We mention in passing that due to the strong damping, weak- 
coupling Redfield-type approaches [11] are not suitable for a description of 
condensed-phase ET processes. 

For very large uj c /A and sufficiently low temperatures, the scaling limit of 
the SB model is reached, where powerful analytical methods are available 
[5,7]. Unfortunately, the scaling regime is of little relevance to ET reactions, 
but in certain parameter regions, progress is nevertheless possible. A famous 
and very successful description of ET reactions concerns the high-temperature 
limit, where the bath behaves classically, leading to Marcus theory [1]. Indeed, 
it can be shown that the classical limit of the SB model directly gives Marcus 
theory [12,13]. Another example is given by the golden rule formula for the 
nonadiabatic ET rate [14,15], which is accurate for arbitrary temperature as 
long as A/a> c <§C 1. In the opposite adiabatic limit of a very slow (classical) 
bath, analytical progress is also possible [16]. A particular focus of our study 
lies on the crossover regime between adiabatic and nonadiabatic ET, which 
represents the most difficult and least understood regime from a theory point 
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of view, in particular at low temperatures. 



In our previous work [8], exactly this regime was studied for symmetric sys- 
tems. Applying the path-integral Monte Carlo (PIMC) scheme of Ref. [8] to 
e 7^ 0, we are able to directly monitor the ET dynamics without any ap- 
proximation. This should be contrasted to alternative numerical routes to this 
problem. For instance, in mixed quantum/ (semi) classical simulations for the 
SB model [17,18,19,20,21], it is hard to justify some of the approximations, and 
hence the accuracy of the results has to be checked case by case. Similar argu- 
ments apply to basis set calculations [22] , and to memory-equation approaches 
[23,24]. Other methods are known to be accurate, but only in either the scaling 
regime or the weak-coupling regime [25,26,27,28]. We conclude that the PIMC 
method represents an excellent computer simulation method for treating the 
crossover regime mentioned above. PIMC is free from approximations, but has 
to deal with the well known dynamical sign problem, because a calculation of 
ET rates requires dynamical information incorporating interfering real-time 
trajectories. Some workers have attempted to use analytical continuation of 
imaginary-time PIMC data [29,30,31], but this process is mathematically ill- 
defined and troublesome to carry out. Here we instead directly proceed in real 
time, where reliable error estimates can be obtained and true nonequilibrium 
preparations can be investigated [8,32]. The outline of the paper is as follows. 
In Sec. 2, we summarize the model and the simulation technique. Results for 
asymmetric ET reactions are presented in Sec. 3. We close with a concluding 
discussion in Sec. 4. 



2 Spin-boson description of ET reactions 



The spin-boson model [5,7] is defined by a coupled system-bath Hamiltonian 
H = H + Hj + Hb- The "spin" corresponding to the TLS is parametrized 
by Pauli matrices a"j, where the |+) (|— )) eigenstate of a z refers to the donor 
(acceptor). In the absence of the solvent, this leads to 

frA he 

Ho = -—° x + y^- (2) 



The environmental modes are modeled by an infinite collection of harmonic 
oscillators, H B) which bilinearly couple to the position of the electron {Hi). 
The solvent-TLS coupling is completely encoded in a spectral density J(co), 
which effectively becomes a continuous function of oo for condensed-phase en- 
vironments, and determines all bath correlation functions that are relevant 
for ET dynamics. Here we take the Ohmic form (1) as a prototype spectral 
density. The bath autocorrelation function for complex time z = t — ir is (with 
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P = 1/ksT) 



L[Z) = J — J{UJ) sinh(^/2) • (3) 

In the classical limit, the reorganization energy HA, which is an integral quan- 
tity describing the overall coupling strength, is the only important bath quan- 
tity. Details of the frequency dependence of J{uS) become crucial only at low 
temperatures. 

Here we study two dynamical properties characterizing ET kinetics. First, 
P(t) = (a z (t)) = (e iHt / h a z e~ iHt / h ) (4) 

gives the difference in occupation probabilities of the donor and the acceptor 
state, with the electron initially fixed on the donor. This quantity then directly 
probes ET dynamics after a nonequilibrium initial preparation, e.g. following 
photoexcitation of an electron into the donor state. Second, the correlation 
function 

C(t) = (a,(0)a,(t))f> = z ~ h * e iHt/h e"^} (5) 

with Z = tr{exp(— PH)} probes the dynamics under an equilibrium prepara- 
tion. At high temperatures, preparation effects are not expected to dramat- 
ically affect the ET dynamics, and therefore a well-defined thermal rate k t h 
should exist [33,34]. This rate then supposedly governs the exponential relax- 
ation of both P(t) and C(t) on a timescale of the order of r re i ax k^. In 
principle, at low temperatures, this assertion can be violated, and then the 
relevant quantity of interest depends on the experimental setup under consid- 
eration. 

In case a unique rate constant exists, it is composed of the forward rate kf 
and the backward rate k^, referring to directed transfers between the donor 
and acceptor states according to 

dP + = {-k f P + + k b PJ)dt , dP_ = (-k b P_ + k f P+)dt 

for times beyond some transient timescale, t > r trans ~ l/u c . Under such a 
rate process, ET dynamics obeys 

P(t > r trans ) = P trans exp(-fctht) + Poo, (6) 
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where = {<J z )p denotes the electronic equilibrium population, 

-Ptrans = [-P('Ttrans) -Poo] 6Xp (&th ^trans ) i 



and A; t h is the total transfer rate kf + k\,. Accordingly, P(i) follows a simple 
exponential decay for times t > r trans . 

Another way to obtain the thermal transfer rate /c t h is via a time-dependent 
equilibrium rate function [33,34]. Assuming that kf and kb are connected via 
a standard detailed balance relation, 

kh = ktexp(—hj3e) , (7) 



for the SB model this function reads [8] 

kit) = ^^ImC(0, (8) 



where we exploited that Eq. (7) corresponds to the assumption of a Boltzmann 
distribution for the equilibrium occupation probabilities, P m = — tanh(^/9e/2). 
This is expected to always hold under the strong damping conditions char- 
acteristic for ET reactions. Indeed, our numerical results below verify this 
assumption to high precision. 

According to the reasoning in Ref. [33], k th is then obtained as the plateau 
value of k{t) if such a plateau exists. However, a careful re-examination of the 
argument in Ref. [33], see also Ref. [35], reveals that this procedure is just 
the limiting case of a more general situation, where kit) decays exponentially 
according to 

k(t > r trans ) = k th e- k ^. (9) 



Therefore, even if Eq. (8) exhibits no clear plateau, a reaction rate can still be 
extracted if the "exponential fit" to Eq. (9) is possible. A concrete example is 
shown in Fig. 1, where the PIMC data for P(t) and k(t) are both consistent 
with the same thermal rate constant /c t h = 0.097A. The exponential fit to 
Eq. (9) is generally more accurate than extracting k^ directly from P(t), where 
transient decays for t < r trans complicate the analysis. The latter approach was 
taken in Ref. [8] whenever no plateau was found for k(t), but then the rates 
could have been obtained to better precision using the exponential fitting 
procedure for k(t). We mention in passing that regarding Fig. 7 of Ref. [8], 
ET rates for T > 2^A//c B indeed closely follow the classical Marcus rate, see 
Eq. (10) below, if one takes a renormalized temperature V = 0.87T (0.78T) 
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Fig. 1. ET dynamics for A/A = 15, e = 0, A/uj c = 0.067, and T = hA/k B . Shown 
are PIMC data for P(t) (top) and k(t) (bottom). Dashed curves are exp(— k t ht) 
(top) and kth exp(— ktht) (bottom), with k t h = 0.097A. Error bars (one standard 
deviation) are smaller than symbol size. 

for A/A = 10 and A/cu c = 1 (A/lu c = 2). Obviously, fits to Eq. (9) can fail, in 
particular when (vibrational) coherence survives, or when r trans and r rc i ax are 
not sufficiently well separated. Such problems will then indicate the invalidity 
of a simple rate picture in these cases. 

Before describing the PIMC method employed here, we briefly review classical 
Marcus theory [1]. Writing the thermal transfer rate in terms of the forward 
rate, k cl = kf(l + e~ ne / kBT ), Marcus theory predicts 



A2 ,g 
f 4 + 7rA 2 /(Aw r ) V Ak B T 



-F*(e)/k B T 



(10) 



with the activation free energy barrier (Marcus parabola) 
F*(e) = ft(e-A) 2 /4A, 



11^ 



and a solvent frequency scale ou T . The u T dependent prefactor in Eq. (10) 
was introduced by Zusman [2] to account for recrossing events. The scale uj t 
seems to follow a power law at low temperatures [8], but in the true classical 
(high-temperature) regime, it is given by uo T = uo c /2 [8,12]. Marcus theory 
(supplemented by dynamical recrossing effects) yields a classical rate constant 
covering the full crossover from nonadiabatic to adiabatic ET. One of its most 
striking features is a non-monotonic behavior of the forward rate as a function 
of asymmetry e. According to Eq. (11), kf exhibits its single maximum in the 
activationless case e = A, but decreases again in the inverted regime e > A. 
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While experiments qualitatively support this picture [36], quantum effects 
(nuclear tunneling) not captured by Marcus theory should play a major role 
in the inverted regime [1]. Such effects are of course included by the PIMC 
simulations below. 

Next we briefly describe our computational technique used to calculate ET rate 
constants (if they exist) under a spin-boson description. A detailed exposition 
of our method can be found in Ref . [8] , and here we only sketch it to keep the 
paper self-consistent. After tracing out the Gaussian bath degrees of freedom 
in Eqs. (4) and (5), one arrives at path-integral expressions for the dynamical 
quantities of interest, where only the TLS degree of freedom is kept explicitly. 
For instance, 

C(t) = Z- 1 J Va a(0)a(t) exp {^S [a} - $[<r]} , (12) 

where the path integration runs over paths <r(z) for the discrete spin variable 
a = ±1 corresponding to the eigenvalues of a z , with the complex time z 
following the Kadanoff-Baym contour ^ = {z: O^t^O^ —ih(3}. S [a] 
denotes the free TLS action, while the effective action due to the traced-out 
bath is captured by the influence functional [5] 

$[cr} = ^J dz J dz'a(z)L(z - z')a(z'), (13) 

7 z'<z 

with integrations ordered along the contour 7. For a numerical evaluation, time 
is discretized, and, moreover, the real-time spins cr(t') and u'{t') on the forward 
and backward branch of the Kadanoff-Baym contour 7 are combined to form 
quantum fluctuations (<r — a') and quasiclassical variables (a + a') which are 
sampled stochastically in the course of the PIMC simulation. Technical details 
about our implementation of this program can be found in Ref. [8]. 

While similar expressions can be derived for the occupation probability, P(t) 
can also be obtained from the same MC trajectory as C(t) (and therefore 
k(t)) by including a correction factor that accommodates the changes in the 
corresponding MC weight [8]. In principle, this allows to simultaneously com- 
pute k(t) and P(t) without significant increase in computing time. However, 
for a strongly biased system, e/A ^> 1, this approach has a rather poor per- 
formance since the respective MC weights refer to almost orthogonal initial 
electronic states. It would then be necessary to run separate simulations to 
extract P(t) and k(t). More serious complications arise from the fact that k{t) 
in Eq. (8) is composed of two competing contributions, namely [1 + cosh(/i / 9e)] 
and lmC(t)/h(3. Since the former increases exponentially with |e|, the latter 
should be very tiny in order to give a finite result. Unfortunately, our simu- 
lations show that the stochastic PIMC error in ImC(t) is rather insensitive 
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Fig. 2. Equilibrium expectation value ((J z )p as a function of e for T = HA/k^ 
(circles) and T = 10ft A/&b (diamonds). The dashed and dotted lines refer to 
— tanh(ft/3e/2). Solvent parameters are A/lo c = 1 and A = 10A. 

to a change of parameters, as long as the time scale of the simulation is kept 
fixed. Hence relative errors become arbitrarily large as ImC(t) — > 0, which ren- 
ders data for k t h numerically unstable. Since the relevant parameter is \h/3e\, 
extraction of k t h from k(t) fails already at rather small |e| for low tempera- 
ture. These numerical problems are unrelated to the dynamical sign problem, 
which did not impose serious limitations in the parameter regime investigated 
here. Therefore, there was no need to employ the powerful but complicated 
"multilevel blocking" algorithm [8] to relieve the sign problem. 



3 Numerical results for asymmetric ET 

Being mainly interested in the effect of the asymmetry e on the ET dynam- 
ics, we have calculated both the occupation probability P(t) and the time- 
dependent rate function k(t) for various values of e at a high and a low tem- 
perature. Specifically, we consider T = 10ftA//cB and T = ftA//c B , and take 
bath parameters where Marcus theory has been shown to be accurate in the 
symmetric case [8], namely A/u c = 1 and A = 10A. This parameter choice 
represents the most interesting crossover region between nonadiabatic and 
adiabatic ET. 

For these temperatures, from closer inspection of Eq. (10), the rate fcth(e) is 
expected to exhibit a single maximum at e = for the higher temperature, 
but two symmetric maxima at e = ±e max with e max < A for the lower one. 



From Eq. (10), it is easy to see that e max is a solution of the equation 
(1 - e/A)[l + exp(ft/3e)] = 2. 

Furthermore, due to Eq. (7), the thermal transfer rate is a symmetric function 
of e. Quantum corrections to the classical rate (10) are then of primary inter- 
est to us, in particular in the inverted regime, |e| > A, for both signs of the 
asymmetry. Before turning to the thermal transfer rate, we show the equilib- 
rium expectation value of a z as a function of e, see Fig. 2. For the parameters 
investigated, we find to high precision (<J Z )[3 = — tanh(/i/3e/2), as has already 
been mentioned above. 

Starting with the high temperature T = lOhA/k-Q, we then calculated the ET 
dynamics for \e\/A < 1.8 while keeping all other parameters fixed. Then the 
difficulties mentioned above did not pose serious limitations. The correspond- 
ing data for /c th , which always follow from exponential fits to k(t), are displayed 
in Fig. 3. Similar to the crossover regime for the symmetric system [8], the 
classical rate (10) again nicely reproduces the exact results. Our PIMC data 
slightly exceed the classical prediction, especially for small |e|. However, these 
deviations can be astonishingly well compensated for by the above-mentioned 
renormalization of the temperature. While we do not have a good theoretical 
argument for the validity of the temperature renormalization, we wish to stress 
that a renormalization of the reorganization energy and/or the activation en- 
ergy is less apt to fit our data over the full asymmetry range, and moreover, 
completely fails to capture the mentioned deviations in the symmetric case. 
Notably, in the high-temperature regime, no significant trace of nuclear tun- 
neling could be observed in the inverted regime. Outside the transient motion, 
ET dynamics exhibits the exponential decay (6), with slower (faster) decay 
than during the transient stage for e < (e > 0), see Fig. 4. 

For the lower temperature T = hA/k-Q, the numerical limitations discussed 
above became quite severe. Extraction of the transfer rate from k(t) could 
only be achieved for |e|/A < 0.6. Nevertheless, for e < 0, k t h could still be 
obtained from Pit), although with poorer accuracy since two fit parameters 
(/c t h and Ptrans) are required. This extends the range of our data to —1.4 < 
e/A < 0.6, and hence includes one of the rate maxima. Again the classical 
rate yields a remarkably accurate prediction for A; t h, see Fig. 3, while the 
temperature renormalization is seen to be limited to the high-temperature 
regime [8]. One might therefore be tempted to conclude that there are no 
significant quantum effects even at low temperatures. However, for e/A < 
—0.2, the ET dynamics shows a qualitatively different behavior. Due to the 
strong equilibrium localization of the electron to the donor state, the fast low- 
temperature initial transient kinetics is able to carry P(r trans ) to a value below 
Poo, causing a subsequent increase of P(t) towards Px,, see Fig 4. Nevertheless, 
the subsequent time evolution is still correctly captured by Eq. (6), albeit the 
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Fig. 3. k th (e) for T = 10hA/k B (top) and T = hA/k B (bottom). The dotted 
(dashed) curve refers to Eq. (10) without (with) a renormalization of the tempera- 
ture, T' = 0.87T. The larger error bars for T = HA/kB and e/A < —0.8 are caused 
by rate extraction from P(t) rather than k(t). 



exponentially decaying excess population with respect to its equilibrium value 
is now negative. This effect nicely illustrates that even outside the transient 
kinetics, only calculating the transfer rate is sometimes not sufficient to fully 
characterize the ET dynamics. 
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Fig. 4. P{t) for T = 10^A//c B (top) and T = hA/k B (bottom). The dot-dashed 
curves denote fits to Eq. (6) with k t h as in Fig. 3 and (top left to top right) 
^trans = 0.41, 1.29, 1.72, and (bottom left to bottom right) -0.11,-0.24,-0.04, 
respectively. For the sake of clarity, only part of the data is shown. 

4 Discussion 

Using the SB model as a description for ET processes, we have calculated 
the thermal transfer rate k t h for a bath setup in the crossover regime, where 
Marcus theory is known to provide a valid description for the symmetric elec- 
tronic system. Focusing on the influence of the energy bias e between the 
redox states, we investigated the electronic dynamics for two different tem- 
peratures, T = hA/kB and T = lOHA/k-B and various values of e up to the 
inverted regime using real-time PIMC simulations. We extended our previ- 
ous approach to extract the thermal transfer rate from the time-dependent 
function k(t), see Eq. (8), which yields a powerful tool both to obtain kth 
and to decide whether a rate description is appropriate in the first place. In 
the absence of electronic or vibrational coherence, this rate also describes the 
exponential decay of the electronic population P(t), provided the timescales 
for transient dynamics r trans and thermal relaxation r rclax are sufficiently well 
separated. 

Our results for the asymmetric electronic system support the rate picture at all 
values of the electronic bias investigated, which extend well into the inverted 
regime. This also applies for the cases where the transient motion shifts the 
electronic population P(t) below its equilibrium value, i.e. at e < and suf- 
ficiently low temperatures, where it is the subsequent deviation of P(t) from 
this equilibrium value which decays exponentially. Furthermore, the validity 
of the classical rate (10) was tested as a function of the asymmetry. Devia- 
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tions from the classical rate expression indicating contributions from nuclear 
tunneling were largely absent, suggesting that they should only play a major 
role at either considerably lower temperatures than T = HA/kB, or for very 
large asymmetry \e\/A. However, we could observe an overall enhancement of 
the rate which can be nicely captured by evaluating the classical rate for a 
renormalized temperature T' < T. This phenomenon was also confirmed for 
the symmetric system, seems to be confined to the high-temperature regime, 
and becomes more pronounced towards the adiabatic limit. Furthermore, at 
low temperatures we found a rapid transient dynamics, followed by a back- 
flow to the donor. This effect shows that a naive rate picture breaks down 
in the low-temperature regime, although the timescale for relaxation is still 
well predicted by the classical prediction. These findings illustrate that despite 
the good accuracy of classical Marcus theory, which is confirmed once again 
by our paper, a fully quantum-mechanical description of the bath's influence 
on the ET process can be important. For the lower temperature, our method 
eventually failed to produce reliable data with increasing |e|, especially for 
e > 0. While there is no obvious cure for the computation of k(t), for P(t) 
these problems seem to mainly stem from our use of a particular PIMC scheme 
optimized for computation of k(t). A direct calculation of P(t) will be free of 
these problems and could therefore allow to significantly extend the range of 
asymmetries. 

This work has been supported by the Volkswagen-Stiftung and by the Deutsche 
Forschungsgemeinschaft . 
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